% computes the theoretical calbiration
% functions of GMS05				cmanoj-16.10.2001
% added the helper functions to
% the main programme 23.4.3.-> Fe.m,FmMFS05.m,fhbp1.m,FN.m,
% sig.m,flp2.m,fhp.m,flp3.m,pv.m,NrFct.m,Fcha.m,
% flp4.m,Ce.m,Ch.m
% Ref: ProcMT manual

% bug fixing details 
% 19.10.2001 in FN.m & FmMSF05.m - corrected
% function pv found to give errors - rectified
% Caution band 5 calibration  still may have problems (fix the sample rate)
% latest date 24.4.2003


function[cplx,frq1] = calib05(sensor,nband,block,probedist,notch,hpb1,gain),
notch = notch';

% -----------declarations-----------

SpRate  = [8.0 4.0 2.0 1.0 0.5 0.25];
blksize = [128 256 512 1024 2048 4096];
DeltaFreq05 = [32.0 1.0 1.0/32.0 1.0/1024.0 (1.0/32768.0)/0.25]; % copied from pmtconst.def
% increment frequency 
% Note for Band5 data, the Deltafreq05 may need to be changed
MaxChnl = 8;

j = sqrt(-1);

PAI = 3.1415926535897;

NFact = NrFct(nband,sensor,probedist,gain); % Normation factor for time series
SampleRate = SpRate(find(blksize==block))*DeltaFreq05(nband);

if nband == 1,
hpb1 = 1;
else,
hpb1 = 0;
end;

hpb1 = 0; % forces high pass filter off

fr = SampleRate;

for spec = 1:block/2,
% 27.07.2002 pvmod is changed ver. of pv 
dphi = pv(nband,fr,sensor) * PAI/180;	% phase correction for ADC
%dphi = pvmod(nband,fr,sensor) * PAI/180;	% phase correction for ADC
cPhase = cos(dphi) + j*sin(dphi);
if sensor <= 2,
cplx(spec) = (Ce(fr,nband,notch,hpb1)*NFact)/cPhase; % Electric field sensors
elseif sensor > 2,
cplx(spec) = (Ch(fr,nband,notch,hpb1)*NFact)/cPhase; % Magnetic field sensors
end;

frq1(spec) = fr;
fr = fr + SampleRate;
end;


%-----------Helper functions
%--------------------------------------------------------------------------
%     Computes the norm-factor for the time series data
%     dynamic cnst of ADC, distances of E-Channels and pregain
%    --------------------------------------------------------------------------
% cmanoj-16.10.2001
% usage NrFct(Band,Sensor,ProDist,PrGn),

% removed TFactor !!
%latest date 5.10.2002


function[factor] =  NrFct(Band,Sensor,ProDist,PrGn),

cDynFactADC = 20.0 / 4194304.0;  
factor = cDynFactADC;     
cFir = 1.25;

if Band ~= 1,
 factor = factor / cFir;
end;

if Sensor <= 2,
factor = factor * 1.0E6 / ProDist;
end;

factor = factor / PrGn;
%--------------------- pv
% returns the phase delay at ADC of GMS05
% 4 channels 4 band -cmanoj 16.11.2001
% refer procMT manual for details
% rev 5/10/2 - poffset for bd 5  set to 0
% latest date 5.10.2002
% Note sensor and Channel number are locked
% 1 to 2 always E channels . 3 to 5 always magnetic channels


function [deltaP] =  pv(band,hz,Sensor),
MaxChnl = 8;
FLP     = [8192 256 8 0.25 0.25/32.0]; % highest frequency in a band
Poffset = [-10 -2.5 -2.5 -2.5 0];  % phase offset for different channels
deltaP = Poffset(band) * (hz/FLP(band)) * (MaxChnl-(Sensor));

%------------------------------------
% function Ce just e field trnfns
% latest date cmanoj-17.10.2001
% latest date 01.11.2002
% Usage [z] = Ce(hz,band,notch,HPB1);
% Ref: ProcMT manual



function[z] = Ce(hz,band,notch,HPB1),

z = Fe(hz)*FCha(hz,band,notch,HPB1);
%-----------------------------------
% function CH just multi plications
% latest date 17.10.2001

function[z] = Ch(hz,band,notch,HPB1),

z = FmMSF05(hz)*FCha(hz,band,notch,HPB1);
%-----------------------------------------
% Fcha just multiplies 
% trfns of different filters
% cmanoj 11.10.2001
% latest date 17.10.2001
%-----------------------------------------
function[z] = FCha(hz,band,notch,HPB1),

z = FN(hz,notch)*flp2(hz,band)*fhp(hz,band)*flp3(hz)*flp4(hz,band);

if HPB1 == 1,
z = z * fhpb1(hz);
end;
%-------------------------------------------
% model function for Electrical
% channel MFS05-EP05 cmanoj
% latest date 17.10.2001

function[z] =  Fe(hz),
GE = 30.0;
j = sqrt(-1);
f1 = 8192.0;     

z1   = 1/GE + j* ((hz/f1)/GE); 
dphi = 15.5 * (hz/f1) * (pi/180);
z2 = cos(dphi) + j*sin(dphi);
z = z1*z2;
%------------------------------------------------
% Phase delay & LP filter correction
% for MFS05 (cmanoj)
% to be used in clibration
% latest date 17.10.2001
% bug z1 fixed 19.10.2001

function[z] =  FmMSF05(hz),
j = sqrt(-1);

f1 = 8192.0;     
PAI = 3.1415926535897;

z1 = 1 - j*1/(hz/4);
z2 = 1 + j*hz/f1;
dphi = 36.42 * (hz/8192)* (PAI/180);
z3 = cos(dphi) + j*sin(dphi);
z = z1*z2/0.8;  
%---------------------------------------------
%notch filter tr. fn. cmanoj 10.10.2001
%latest date 17.10.2001
%K = Notchi(find(IsNotch>0)) use this for sending the active notches

function[z] = FN(hz,notch),
j = sqrt(-1);

aaa = length(notch);

RQ = 1/5.0; % standard for GMS05
b = 0;
for i = 1:aaa,
b = b + sig(hz,notch(i));
end;
b = b * RQ;
z = 1 + j*b;
%-----------------------------------------------
% highpass filter 
% % latest date 24.03.2003
% j=sqrt(-1) added 24.3.2003

function[z] = fhp(hz,band),
j=sqrt(-1);

HPFP = [0.25 0.25 0.25 0.00048828125 0.00048828125];

z1 = 0 + j* (hz/HPFP(band));
z2 = 1.0 + j* (hz/HPFP(band));
z = z2/z1;
%--------------------------------------------------
% high pass filter for band 1
% latest date 17.10.2001

function [z] = fhpb1(hz),
j = sqrt(-1);

f = hz;

z1 = 1 + j* (-1/(f/800));
z2 = (1.0 - 1.0/((f/800)*(f/800))) + j * (-1/(f/800));
z = z1*z2;
%---------------------------------------------------
% Low pass filter for band  2 & 4
% latest date 17.10.2001

function[z] = flp2(hz,band),
FLP = [8192 256 8 0.25 0.25]; 	%for bd # 1 2 3 4 (band 5 ?)
j = sqrt(-1);
z = 1 + j* (hz /FLP(band));
%-----------------------------------------------------
% low pass filter for band 3
% latest date 17.10.2001

function[z] = flp3(hz),

f1 = 8192.0;
j = sqrt(-1);
a  = 0.63333;

z = (1-(hz/f1)*(hz/f1)) + j*a*(hz/f1);
%-------------------------------------------------------
% Low pass filter for 2 & 4
% latest date 17.10.2001

function[z] = flp4(hz,band),
FLP = [8192 256 8 0.25 0.25]; 	%for bd # 1 2 3 4
j = sqrt(-1);
z = 1 + j* (hz /FLP(band));
%-------------------------------------------------
% divn 
% used with calibration of GMS05 channels
% latest date 17.10.2001

function[divn] = sig(hz,hzn),

GRTNUM = 999.99;

RealTmp = (hz/hzn)*(hz/hzn);
if RealTmp ~= 1.0,
divn = (hz/hzn) / ( 1.0 - RealTmp);
else
divn = GRTNUM;
end;
%-------------------------------------------------